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Abstract. 

We describe a computational framework linking Uncertainty Quantification (UQ) methods for 
continuum problems depending on random parameters with Equation-Free (EF) methods for per- 
forming continuum deterministic numerics by acting directly on atomistic/stochastic simulators. 
Our illustrative example is a heterogeneous catalytic reaction mechanism with an uncertain atom- 
istic kinetic parameter; the "inner" dynamic simulator of choice is a Gillespie Stochastic Simulation 
Algorithm (SSA). We demonstrate UQ computations at the coarse-grained level in a nonintrusive 
way, through the design of brief, appropriately initialized computational experiments with the SSA 
code. The system is thus observed at three levels: (a) a fine scale for each stochastic simulation 
at each value of the uncertain parameter; (b) an intermediate coarse-grained state for the expected 
behavior of the SSA at each value of the uncertain parameter; and (c) the desired fully coarse-grained 
level: distributions of the coarse-grained behavior over the range of uncertain parameter values. The 
latter are computed in the form of generalized Polynomial Chaos (gPC) coefficients in terms of the 
random parameter. Coarse projective integration and coarse fixed point computation are employed 
to accelerate the computational evolution of these desired observables, to converge on random sta- 
ble/unstable steady states, and to perform parametric studies with respect to other (nonrandom) 
system parameters. 

Key words. Uncertainty quantification, Equation-Free, generalized Polynomial Chaos, stochas- 
tic Gillespie algorithm, multiscale 

1. Introduction. The temporal evolution of many engineering systems can be 
described through continuum models (typically Ordinary or Partial Differential Equa- 
tions, ODEs or PDEs) depending on random parameters and/or initial/boundary 
conditions. A number of methods have been proposed to study the evolution of the 
probability distribution of the solutions of such random problems (the so-called un- 
certainty quantification or UQ). Among early exploration approaches in this area are 
Monte-Carlo based methods ^JE^EIE], which normally require a large number of 
ensemble realizations to achieve convergence. For slightly perturbed systems, whose 
random parameters can be described as small fluctuations around average values, 
perturbation methods may be utilized 000. Such methods cannot, however, be 
used to treat more generic random systems, exhibiting large parametric uncertain- 
ties, and even if used they can only obtain information on low-order statistics. To 
overcome this difficulty, moment-closure techniques (e.g., were tried to study dy- 
namics of statistical moments of solutions for which the small-uncertainty assumption 
in the perturbation method may be relaxed. A significant difficulty with this aprroach 
lies in deriving an accurate closure for high-order statistical moments of the solution 
distributions; this is extremely difficult, especially for nonlinear systems. 

In recent years an alternative approach for UQ, the stochastic Galerkin method, 
has received considerable attention as an approach to solving ODEs or PDEs with 
random parameters. This approach originated from the work of Wiener |5] in con- 
structing multiple stochastic integrals (also known as Homogeneous Chaos) to repre- 
sent functionals of Wiener processes. This idea was then utilized in JU] to express 
solutions of random systems in terms of Hermitc polynomials of Gaussian random 
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variables (i.e., Wicner-Hermite Polynomial Chaos). Projections of these solutions on 
the Hcrmitc polynomial basis are deterministic, and can be readily numerically solved 
after truncation through a Galcrkin method. Early works in applying the method in 
engineering systems include |10j . where random static and dynamical problems in 
structural analysis were investigated. The method was subsequently applied in uncer- 
tainty quantification of various physical systems including (but certainly not restricted 
to) porous media fluid dynamics ^2] and chemical reactions [Ej. The applica- 
tion of the Wiener-Hermite Polynomial Chaos was recently extended in |141 1151 
to a general situation, where the Askey scheme was used to construct generalized 
Polynomial Chaos (gPC) with respect to other continuous and discrete probability 
measures. Other developments in the application of the method include pieccwisc 
function representation |17j , non-orthogonal expansions |18j and Wiener-Haar wavelet 
expansions ^U]- O ne advantage of the stochastic Galerkin method relative to Monte 
Carlo simulation is that it can reduce a random system to a deterministic one, often 
with considerably fewer degrees of freedom. Moreover, this method - along with its 
extensions and ramifications - can handle more general uncertainty problems such as 
large fluctuations, multimodal and discontinuous probability distributions, etc. 

In order to implement the stochastic Galcrkin method, we must be able to derive 
differential equations for the projection of solution distributions onto gPC bases (i.e., 
ODEs for the evolution of gPC coefficients). When explicit equations for these gPC 
coefficients are not available in closed form, a unified approach, combining equation- 
free techniques |20l 1211 |2~2] and the stochastic Galerkin method in a nonintrusive 
way has been proposed. This is equation-free uncertainty quantification |23j . which 
can evolve gPC coefficients in time, and has been used to perform steady state and 
limit cycle bifurcation analysis of the gPC coefficient equations without needing these 
equations in closed form. This approach is based on short bursts of direct simula- 
tion (ensembles of such simulations distributed over the random parameter(s)) as an 
inner, "fine scale" simulator; these bursts are used to numerically estimate the neces- 
sary information at the "coarser" gPC coefficient level (e.g., temporal derivatives of 
gPC coefficients) on demand. The main assumption underlying this approach is that 
the long-term dynamics in gPC coefficient space lie on a low-dimensional, attracting 
manifold, which can be parametrized by only a few leading-order gPC cocffiiccnts 
in the appropriate orthogonal polynomial chaos basis (see |14)V In this way, model 
reduction can be achieved by using gPC expansions. 

When we study continuum models of chemically reacting systems (e.g. ODEs 
for the evolution of reactant concentrations in a stirred tank reactor, or for the evo- 
lution of coverages on catalytic surfaces), UQ techniques can be used to study the 
effect of uncertainty in kinetic or operating parameters on the overall reactor be- 
havior (e.g. steady state concentrations, reaction rates etc.). In many problems of 
contemporary interest, however, such continuum differential equations are not avail- 
able in closed form; instead, we are given a description of the chemical kinetics at 
an atomistic/stochastic level and the uncertainty enters now in the parameters of 
the atomistic/stochastic simulation itself (e.g. in certain transition probabilities). 
The normal steps for uncertainty quantification modeling would involve (a) deriva- 
tion of closed continuum kinetic equations for the coarse-grained observables of the 
atomistic simulation (averaged concentrations, averaged coverages)- these equations 
need to explicitly express coarse-grained parameter uncertainty in terms of the fine 
scale, atomistic level parameter uncertainty; and (b) the application of traditional, 
continuum UQ techniques on these continuum kinetic equations. Here we will show 

2 



how to circumvent this "intermediate" derivation of coarse-grained, continuum kinetic 
equations. 

We will implement a computational framework that performs UQ computations 
for these unavailable, random continuum kinetic equations by acting on the fine scale, 
atomistic/stochastic simulator directly. This uncertainty quantification procedure for 
stochastic chemical reaction models with uncertain parameters involves three distinct 
system levels: a fine scale, microscopic simulator for fixed values of the uncertain 
parameter (s); an intermediate coarse-grained scale, where continuum kinetic ODEs 
exist in principle for certain system observables, still for fixed values of the uncer- 
tain parameter(s); and, a final, fully coarse-grained scale of gPC coefficients for the 
distributions of the solutions of these (unavailable) closed form continuum equations. 
The equation-free machinery is thus used at two successive levels: first to circumvent 
the derivation of closed form continuum kinetic equations; and then to circumvent 
the derivation of explicit equations for the evolution of the gPC coefficients for the 
solutions of these equations with uncertain parameters through a Galerkin procedure. 

The paper is organized as follows: We start in Section [21 briefly recalling the 
formulation of the stochastic Galerkin method. Equation-free techniques are com- 
bined with this stochastic Galerkin method in Section |21 to study random systems 
without explicitly available governing equations. The approach is illustrated through 
the standard model for the A + 1/2B-2 — > AB reaction, which has been used as a 
simplified description for the catalytic oxidation of CO [21] ■ The fine-scale simulator 
is chosen to be the Gillespie Stochastic Simulation Algorithm (SSA). Computational 
results are presented in Section^where we also demonstrate the efficiency of using the 
Gauss-Lcgcndrc quadrature for computing gPC projections. Random stable/unstable 
steady-state solutions are also computed in Section^] and continuation algorithms are 
implemented to explore the effect of variation of a system (nonrandom) parameter on 
the statistics of the random steady state. We conclude with a brief summary and 
discussion in Sectional 

2. The Stochastic Galerkin Method. Consider a system whose evolution in 
time is governed by the differential equation 



the system state is x c (u> c ,t), where lo c is an element in the sampling space Q c . The 
subscript c stands for "coarse" , since we consider this to be the coarse-grained descrip- 
tion of an atomistic simulator (an SSA simulator in our illustrative example below) . 
The solution of the above equation can be described by an expansion in an L 2 space 
with a generalized polynomial chaos basis ^i($(uj c )), i.e., 



(2.1) 



— - = f(x,w c ), x c (uj c ,0) = x cfi (uj c ); 



p 



(2.2) 




The projections x l cc (t) are determined by 



(2.3) 



< x c (u c ,t),Vi(£(Lj c )) > 

<* 1 (€(Wc)),* i (€(Wc))>' 



where the inner product < • , • > is defined as 



(2.4) 
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(where p(£) is a probability measure of £) for any functions q(x(£)) and g(x(£)) in 
the L 2 space. The subscript "cc" is indicative of a second level of coarse-graining: 
these coefficients describe the statistics of the distribution of solutions for the already 
coarse-grained state x c (ui c , t) over the sampling space. 

Through a Galcrkin projection, equations governing gPC coefficients x l cc (t) are 
obtained as 

dx\ c _ < /(E^Qg'ccWttite)).*^) > 

dt <**(€),*<(€)> 

(2.5) i = 0,2,---,P, 

with x l cc (0) = ^^[gj'^Hj^ ■ The above equation can be formally compactly rewritten 
as 

(2-6) ^ =H (X CC ), 
where 

(2.7) X cc = (x cc , x cc , • ■ ■ , x cc ) 
and 

(2.8) H = (ho, hi, ■ ■ ■ , hp) T . 



Here 



hi(X c 



</(Ef=o4c(*)*i(0) J ^(0> 



<MZ),MS)> 

(2.9) i = 0,2,---,P. 

If d * t cc — in the long-time limit, then Eau. (|2.tj|) has a steady state, which can 
be used to obtain the probability distribution of the random steady state of 1(2.1(1 . 
The explicit derivation of Equation 1)2.6(1 is a challenging problem if 1(2.1(1 is a set 
of strongly nonlinear equations; pseudospectral approaches (see e.g. (13p provide a 
possible alternative. 

3. Equation-Free Computation for Random Dynamical Systems with- 
out Explicit Governing Equations. Equation-free methods have been applied in 
recent years to investigate solutions of non-random macroscopic systems whose evo- 
lution equations are not explicitly available [2(JI 1241 HH1 1221 125( . The approach can, 
in principle, provide clear scenarios of the coarse-level evolution and its parametric 
dependence while requiring only short-time bursts of evolution with the micro-level 
simulators; in effect, it is a framework of accelerating the extraction of information 
from the microscale simulation through judicious design of computational experiments 
and processing of their results. 

The equation-free approach utilizes the so-called coarse time-stepper as its basic 
element; this time-stepper consists essentially of three components: lifting, micro- 
simulation, and restriction. Lifting is a protocol that transforms a coarse-level state 
to consistent fine-level states; restriction is the converse of lifting. Note that the lifting 
will, in general, not be a one-to-one transformation, since fine-scale states have far 
more degrees of freedom than their corresponding coarse-grained descriptions; this is 
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a vital step in EF computations, and a good code should test, on line, that different 
realizations of the lifting protocol do not affect the coarse-grained computational 
results. 

This approach was applied recently in conjunction with the stochastic Galcrkin 
method to study random dynamical systems |23| . The fine scale state, in this case, is 
a (large enough) ensemble of system realizations; at each time snapshot, this ensemble 
can be represented by its projection onto a generalized polynomial basis; if a low order 
gPC truncation provides an accurate representation, these first few gPC coefficients 
constitute the coarse-grained description. In principle, there exist differential equa- 
tions for these projections, which can be derived, implemented in a computer code, 
and then solved using numerical methods. However, it may be difficult or even impos- 
sible to derive these equations; equation- free methods were utilized as an alternative 
way of solving them without deriving them first. 

Assuming that the coarse-grained evolution is smooth, one can use techniques like 
projective integration j5H]|57] to accelerate the time-evolution of the gPC coefficients 
using fine scale (ensemble realization) simulations over only relatively short time seg- 
ments. This is accomplished by observing the evolution of the ensemble MC runs on 
their gPC coefficients (through restriction), and use of these data to locally estimate 
the time-derivative of the coarse-grained description (the local time derivative of the 
gPC coefficients). This information is then "passed" to traditional continuum nu- 
merical initial value problem solvers (ranging from the simple explicit forward Euler 
to Runge-Kutta type or even implicit integrators) that "project" the coarse-grained 
state forward in time j^E2]. One thus solves the initial value problem for the coarse- 
grained description with the necessary quantities (the gPC local time derivatives) 
obtained not through a function evaluation from a closed-form model, but through 
processing the results of "judicious" numerical experiments with the fine scale code. 
Beyond coarse projective integration, when Equation l|2.5(l has a stationary state, one 
can turn the coarse time-stepper into a fixed-point operator and use matrix-free meth- 
ods such as Newton-GMRES j^Hl to compute stable/unstable coarse-grained steady 
states. These correspond to random steady states of the original random equations. 

The equation- free technique in |23| . however, requires the evolution equations 
for the random dynamical systems to be explicitly available. In the case that these 
equations are not explicitly available, we now show how to exploit fine-scale mod- 
els underlying these differential equations. We now have two successive lifting (and, 
correspondingly, two successive restriction) levels. To obtain a numerical represen- 
tation for evolution of the desired, "doubly" coarse-grained representation (the gPC 
coefficients), the fine-scale states of these fine-scale models are first restricted to an 
intermediate coarse level: we obtain individual states in the ensemble of random ODE 
solutions. At a second level of restriction, the entire ensemble of these states is used 
to compute their gPC coeffcients. We view the level of random differential equation 
as the "intermediate coarse" scale and the level where the gPC projections reside as 
a desired "fully coarse" scale. 

The interaction between intermediate coarse and the fully coarse scale states is 
embodied in equations (|2.2|) and (|2.3|1 . which correspond to lifting and restriction, 
respectively, between these two scales. The interaction between fine and "intermedi- 
ately coarse" states at fixed values of the random parameters £ is described by two 
new lifting and restriction mappings: fi (lifting) and M. (restriction), 



(3.1) 
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Fig. 3.1. Schematic of coarse projective integration for multiscale systems with uncertainty. 



(3.2) 



x c (£,t) = M(X f (Z,t)). 



Recall that the fine-scale states Xf are characterized by many more degrees of freedom 
than their intermediate-scale x c counterparts. As the fine-scale states are restricted 
to the intermediate-scale level, these additional degrees of freedom are eliminated. 

Coarse projective integration for random systems whose coarse-scale equations 
are not explicitly available can be summarized in the following steps: (see Fig. 13.111 

1. Generate an ensemble of (intermediate) coarse-scale states based on initial 
values of gPC coefficients by 12.211 . 

2. Generate an ensemble of corresponding fine-scale states consistent with each 
element of the ensemble of (intermediate) coarse-scale states by <|3.1[) . In the 
case that (intermediate) coarse-scale states are mean fields, the average of 
these fine-scale states should be equal to the prescribed (intermediate) coarse 
state. 

Evolve fine-scale states via fine-scale (in the example below, SSA) simulators. 
Restrict the fine-scale states to the (intermediate) coarse level and obtain an 
ensemble of these coarse-scale states by (13.211 . 

Further restrict (intermediate) coarse-scale states to the desired fully coarse 
level, to obtain gPC coefficients by l|2.3J) . 

Perform the above two steps successively, and use the results to estimate the 
temporal derivative of the fully coarse observables (gPC coefficients). The 
data collection time is dictated by the separation of fast/slow time scales 
we assume prevails at the coarse grained level; also by the noise of gPC 
coefficients brought in by the SSA at the fine level and the approach used to 
compute these coefficients at the coarse level. 

Project forward in time -using a continuum numerical integrator, such as 
forward Euler- to obtain the fully coarse observables at a later time. Go back 
to Step 1. The selection of the projective time step so as to retain overall 
stability of the projective method is discussed in |2f>| . 
In equation-free computations of random coarse-scale steady states, we use the 
time-stepper for the fully coarse scale states X cc (the gPC coefficients), <&t-, to con- 
struct a fixed point equation 



3. 
4. 



G. 



7. 



(3.3) 



X cc = $> T {X CC ) 



The operator &t involves repeated lifting and repeated restriction procedures across 
two scale gaps, as illustrated in Fig. 13.21 The solution of Equation (|3.3|) can be 
attemped either by direction iteration, by Newton's method with numerically esti- 
mated Jacobians (for a small number of gPC coefficients) or, more systematically, by 
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Fig. 3.2. Fixed-point formulation for coarse random steady states. 



matrix-free fixed point algorithms such as Newton-Krylov GMRES (for large num- 
bers of gPC coefficients) [2HJ - Once the fully coarse-scale steady states are available, 
ensemble realizations and thus probability distributions of coarse-scale steady states 
can be immediately obtained by l|2.2[l . 

4. Numerical Results. Our illustrative example involves the A + 1/2B 2 — > AB 
reaction, which can be used as a caricature description of the CO oxidation on a 
Pt catalyst surface |2U- Our (intermediate) coarse-scale observables are the mean 
coverages of rcactants and the vacant catalyst sites. The fine scale states in each 
detailed simulation consist of the numbers of sites occupied by each reactant and the 
vacant sites. In the limit of very large systems, the (intermediate) coarse scale ODEs 
for mean coverage based on the particular kinetic mechanism would be 

dd A /dt = a>0* - -y9 A - k r e A e B 

(4.i) de B /dt = pel - k r e A e B 

(where 9 A , 9 B and 9* are mean coverages of reactants A, B and vacant sites, respec- 
tively). We will not use these equations in our computations that follow (which are 
performed for finite size systems) . 

At the fine scale description level, there are four elementary reaction steps: 



A{g) 



-62(5) + *i + *j — > + B*j 



A*,i A(g) 



k, 



(4.2) A M + B* tj AB(g 



here (g) refers to gas phase reactants, a vacant site, and A(B)*^j) are the ab- 
sorbed reactants on the surface. At the continuum limit, the rates of four reactions are 
given respectively by r x = a9*N tot , r 2 = /39lN tot , r 3 = ^9 A N tot and r 4 = k r 9 A 9BNt 



tot J 



where N to t is the number of sites on the reacting surface. For our finite-size system, 
the four reaction rates at a given time t (based on which the reaction probabilities are 
computed) would be n = aN*, r 2 = r 3 = lN A and r 4 = jf^N A N B , 

where N A , N B and A^* are numbers of sites taken respectively by A, B and vacant 
slots at time t. 
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Our implementations of this chemical reaction scheme follows the Gillespie al- 
gorithm (291 I5U| . The index of reaction that will take place next depends on the 
random variable p\, uniformly distributed over the domain [0, 1]: reaction j will oc- 
cur if J2i=i r i/ Si=i r i — Pi — Yji=i r i/ Si=i r i- The time at which this reaction 
takes place is —ln(p2)/ Ylj—i r i, where pi is also a uniformly distributed random vari- 
able over [0, 1]. For each realization in the fine- level SSA, a stoichiometric matrix is 
then used to keep track of the changes in the numbers of the reactants and vacant 
sites over time. 

In the simulations that follow, the parameters a, 7, and k r arc considered known, 
and set respectively to 1.6, 0.04 and 4; the uncertain parameter is (3: (3 = 6.0 + 0.25£, 
where £ is a random variable uniformly distributed over [—1,1]. Legendre polynomial 
chaos is chosen as the basis for the random, fully coarse-scale states. The highest 
truncation order for these fully coarse observables is chosen as 3 (P = 3). A total 
ensemble of 40,000 realizations are used in lifting and restriction between the inter- 
mediate coarse (mean coverages) to the fully coarse (gPC coefficients) level. Consis- 
tent with each of these 40,000 ensemble realizations, 1000 fine-scale realizations were 
simulated; each of them uses 200 2 (Nt o t = 200 2 ) sites on the catalyst surface. The 
fine- (intermediate) coarse restriction Ai consists in taking the average coverage over 
these 1000 fine-scale realizations corresponding to each (intermediate) coarse observ- 
able (i.e., mean coverage of each reactant). We perform successive fine- intermediate 
coarse and then intermediate-fully coarse restrictions at 40 successive coarse-scale time 
steps (Ai c = 0.01s). We then use least-squares fitting to estimate temporal deriva- 
tives of the gPC coefficients based on values at the last 5 time steps. These numerical 
derivatives are then used (in a simple, forward Euler projective scheme) to calculate 
gPC states after a relatively large time step (At cc = 0.8s). Figure |4~T1 shows evolution 
of the mean coverage of reactant A as a function of the random variable £ in the 
time domain; the "empty" time intervals in the plot are the projective forward Euler 
"jumps", over which we do not simulate. Figures l4~2l and l4~3l contrast the evolution of 
the fully coarse-scale observables (gPC coefficients), computed by ensemble average 
from direct Monte Carlo simulation of the coarse-scale ODEs (14.1(1 . and also com- 
puted through projective integration of the two-scale-gap system; Fig. 14.41 compares 
the standard deviations of mean coverages computed from the gPC coefficients in the 
two approaches. The results indicate good agreement between projective integration 
computations and true evolutions at the fully coarse gPC level. 

In an attempt to further accelerate the computation, the Gauss-Legendre quadra- 
ture was used to approximate the inner product in the intermediate coarse to fully 
coarse restriction 1(2. 3|l . In the corresponding fully- to- intermediate coarse lifting, only 
coarse states corresponding to values of £ at the Gauss-Legendre points were gener- 
ated. Figures 14.514.71 show coarse projective integration results using this method, 
which utilizes only 200 coarse-scale realizations (rather than 40,000 ones). To display 
the effectiveness of this technique, we implemented the original Monte-Carlo simu- 
lations in the intermediate level by using a total of 200 realizations; the results are 
shown in figures |4~H1 and 1431 With such a small number of realizations in a Monte 
Carlo simulation, while the zeroth-order gPC coefficients of mean coverages of reac- 
tants are well captured, other, higher order coefficients deviate significantly from their 
true trajectories even at the beginning of the projective integration. This is clearly 
due to the lack of sufficient realizations of mean coverage when ensemble averaging 
is used to compute the corresponding gPC coefficients. Lifting only around Gauss- 
Legendre quadrature points can enhance the effectiveness of our approach to simulate 
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Fig. 4.1. Evolution of 6 a as a function of the random variable £ through coarse projective 
integration 
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Fig. 4.2. Evolution of the zeroth-order gPC coefficient of the mean coverages computed by 
coarse projective integration, with Ne=40,000 randomly selected £ values; symbols: CPI results; 
lines: gPC coefficients obtained via Monte Carlo simulation of the coarse ODEs. 

the evolution of gPC coefficients by using a much smaller number of intermediate 
coarse-level realizations than that required by the standard Monte-Carlo simulation. 
Computing the gPC coefficient evolution by sampling only close to qudarature points 
has been proposed (for Gauss-Hermite quadrature, in a problem of approximating PC 
coefficients for random temperature distribution) in 

We also use a matrix-free Newton-Krylov GMRES method to converge on the 
(deterministic) stable/unstable fully coarse steady states of the gPC coefficient de- 
scription, out of which random (intermediate) coarse steady state distributions can 
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Fig. 4.3. Evolution of additional (higher- order) gPC coefficients of the mean coverages com- 
puted by coarse projective integration with Ne=40,000 randomly selected £ values; symbols: CPI 
results; lines: gPC coefficients obtained via Monte Carlo simulation of the coarse ODEs. 
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Fig. 4.4. Evolution of the standard deviations of the mean coverages computed by coarse pro- 
jective integration, with Ne=40,000 randomly selected f values; symbols: CPI results; lines: Monte 
Carlo simulation of the coarse ODEs. 



be obtained by 12. 2 [I . In this computation, j3 =< j3 > (1.0 + 0.05£), where £ is 
again a uniform random variable in [—1,1]. In Fig. 14.101 dashed lines are true en- 
velopes of random coarse-scale steady states (computed by setting /3 = 1.05 < /? > 
and f3 = 0.95 </?>), and solid lines are the derterministic steady states when 
/?=</?>. Error bars and stars represent respectively ranges and means of ran- 
dom steady states of mean coverages computed using the matrix-free, time-stepper 
based Newton-Krylov-GMRES method. Clearly, this equation-free fixed-point com- 
putation can correctly reproduce the random steady states. Again, values of £ at 
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Fig. 4.5. Evolution of the zeroth- order gPC coefficient of the mean coverages computed by 
coarse projective integration for Ne=200 5 values at Gauss- Legendre points; symbols: CPI results; 
lines: gPC coefficients obtained via Monte Carlo simulation of coarse ODEs. 




time 



Fig. 4.6. Evolution of additional (higher- order) gPC coefficients of the mean coverages com- 
puted by coarse projective integration with Ne=200 £ values at Gauss-Legendre points; symbols: CPI 
results; lines: gPC coefficients obtained via Monte Carlo simulation of the coarse ODEs. 



Gauss-Legendre points are used to generate realizations of intcrmediately-coarse-scale 
observables (mean coverages) and implement fixed-point computations (Fig. I4.11|) . 
The random steady states can be accurately captured by this technique as well, while 
the computational load decreases significantly. The approach has been linked with a 
continuation algorithm to trace the bifurcation diagrams of the random steady states 
as a function of < (3 >; observe that unstable random steady states can thus be 
computed, and bifurcation points (such as turning points of random steady states) in 
parameter space can be approximated (see Figures T4 . 1 Ul and 14 . 1 If) . 
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Fig. 4.8. Evolution of the zeroth-order gPC coefficients of the mean coverages computed by 
coarse projective integration for Ne=200 randomly selected £ values; symbols: CPI results; lines: 
gPC coefficients obtained via Monte Carlo simulation of coarse ODEs with Ne=40,000. 



5. Summary. Equation-free methods were used at two successive layers, in this 
paper, to enable uncertainty quantification computations on models of reacting sys- 
tems for which no coarse-grained, continuum description is available in closed form. 
Coarse projective integration was used to accelerate the computation of transient 
dynamics of the problem solution distributions (at gPC coefficients level). Random 
stable/unstable steady state computation, and their parametric/bifurcation analysis 
at this gPC coefficient level was also demonstrated. Gauss quadrature rules were used 
to effectively reduce the computational load while preserving the accuracy of results 
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Fig. 4.9. Evolution of additional, higher-order gPC coefficients of the mean coverages computed 
by coarse projective integration with Ne=200 randomly selected £ values. Red circles and dots 
represent the 3rd- and lst-order gPC coefficients of 8a obtained by CPI, respectively. Jumps during 
the course of their evolution are caused by insufficient number of coarse realizations of coverages 
used to perform ensemble average when computing the gPC coefficiens. Clustered green and blue 
symbols in the lower part of the figure stand for gPC coefficients of reactant 8g and vacancy 8* , 
respectively. Lines represent gPC coefficients obtained via Monte Carlo simulation of the coarse 
ODEs with Ne=40,000. Red, green and blue solid lines stand for respectively the lst-order gPC 
coefficients of 8a, 8b and 8*. Red, green and blue dot-dashed lines stand for respectively the 2nd- 
order gPC coefficients of 8a, 8b and 8* . Red, green and blue dotted lines stand for respectively the 
3rd-order gPC coefficients of 8a, 8b and 8* . 




Fig. 4.10. Bifurcation diagram of the random model steady states with respect to the mean 
< P > of the uncertain parameter distribution. Coarse fixed point computation and continuation 
with Ne=40,000 randomly selected £ values; red and green objects: stable steady states; blue and 
magenta objects: unstable steady states. See text for the inset. 
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FlG. 4.11. Bifurcation diagram of the random model steady states with respect to the mean 
< (3 > of the uncertain parameter distribution. Coarse fixed point computation and continuation 
with Ne=200 £ values at Gauss- Legendre points; red and green objects: stable steady states; blue 
and magenta objects: unstable steady states. See text for the inset. 

of the equation-free results. The method should be applicable in UQ for a wider class 
of reacting problems, beyond those described through a fine-scale SSA simulator, for 
which no closed-form kinetic equations are available. We believe that the approach 
can still serve in cases where we know how to describe uncertainty in the microscopic 
simulation parameters, but we cannot easily translate that in uncertainty descriptions 
for parameters at the coarse-grained level. 
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